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Abstract 

A new class of high-order monotonicity-preserving schemes for the numerical solution of conserva¬ 
tion laws is presented. The interface value in these schemes is obtained by limiting a higher-order 
polynomial reconstruction. The limiting is designed to preserve accuracy near extrema and to work 
well with Runge-Kutta time stepping. Computational efficiency is enhanced by a simple test that 
determines whether the limiting procedure is needed. For linear advection in one dimension, these 
schemes are shown to be monotonicity-preserving and uniformly high-order accurate. Numerical 
experiments for advection as well as the Euler equations also confirm their high accuracy, good 
shock resolution, and computational efficiency. 
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Introduction 


We consider higher-order schemes (at least third-order) for the numerical solution of the Euler 
equations. Typical solutions to these equations have smooth structures interspersed with disconti¬ 
nuities. The challenge is to develop schemes that are highly accurate in smooth regions and have 
sharp nonoscillatory transitions at discontinuities. 

Achieving this dual objective remains a daunting task. Among the first attempts, Colella and 
Woodward [2] introduced a piecewise parabolic method (PPM), which employs a four-point centered 
stencil to define the interface value; this value is then limited to control oscillations. The centered 
stencil, however, results in a scheme with a large dispersion error, and the limiting procedure causes 
accuracy to degenerate to first-order near extrema. 

The essentially nonoscillatory (ENO) schemes of Harten et al. [3] were developed via a different 
line of thought. In these schemes, an adaptive stencil is used to select the “smoothest"' data, 
thereby avoiding interpolations across discontinuities. While an adaptive stencil does avoid spurious 
oscillations near discontinuities, it does not make use of all the available data. The weighted-ENO 
(WENO) schemes by Liu et al. [11] and Jiang and Shu [9] make better use of the available data 
by defining the interface value as a weighted average of the interface values from all stencils. The 
weights are designed so that in smooth regions the scheme nearly recovers a very accurate interface 
value using all stencils but, near discontinuities, it recovers the value from the smoothest stencil. 
The WENO schemes, however, are still diffusive: they smear discontinuities nearly as much as the 
ENO schemes. 

In this paper, we follow the limiting approach. The interface value is defined by a five-point 
stencil. As a result, the leading error is dissipative, and the dispersive error is considerably smaller 
than that of the four-point stencil. This interface value also combines well with the three-stage 
Runge-Kutta time stepping. Similar to PPM, oscillations are controlled by a limiting procedure. 
The key differences, however, are that our limiting is designed to preserve accuracy near extrema 
and to work well with Runge-Kutta time stepping. The resulting scheme is accurate in smooth 
regions, resolves discontinuities with high resolution, and is also efficient. 

Note that a piecewise linear scheme of this type was presented by Huynh [(>]. Extensions to 
piecewise parabolic schemes were presented by Suresh [18] and Huynh [7]. The present scheme 
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incorporates these ideas within a Runge-Kutta time integration framework. 


In §1, the spatial discretization and the Runge-Kutta time integration are reviewed. Section 2 
describes the reconstruction procedure, which is the key feature of our scheme. Extensions of this 
scheme to systems of equations and multi-dimensions are dealt with in §3. Numerical experiments 
appear in §4. Finally, conclusions are presented in §5. 

1. Discretization 

For simplicity, we describe the methods for the advection equation with constant speed a, 


+ 

£ 

II 

o 


u(x y 0) — Uq ( X ) 

(1.16) 


where f is time, x is distance, and uo(x) is the initial condition. For the moment, the solution is 
assumed to be periodic in x so that no boundary conditions are needed. 

Let Xj be the cell center of a uniform mesh, x j+1 / 2 the interface between the j-th and j + 1-th 
cells, and h the cell width. Denote by Uj(t) the cell average of u at time 

1 f X J + 1/2 

= - u(x,t)d.x. (1.2) 

n 1/2 

Integrating (1.1a) over the cell 1 / 2 ] yields 

^ [ w ( ar j + l/2'0 “ m ( 3 j-1/2^)] — (L3) 

At time t n = nr where r is the time step, assume that we know vj which approximates Uj(t n ). 
We wish to calculate . For simplicity of notation, we omit the superscript n when there is no 
confusion, e.g., v 3 denotes vj. 

An approximation to u(xj +1 / 2 , t n ) is called the interface value and is denoted by t^+i/ 2 * The 
calculation of the interface value from the known cell averages is accomplished in two steps. In 
the first or reconstruction step, nonoscillatory approximations of u(#j+i/ 2 ^ n ) 1° the left and right 
sides of the interface x 3 +\/2 denoted by and are constructed. This step determines 

the scheme’s order of accuracy and is the main concern of this paper. In the next or upwind step, 
the interface value is determined by the wind direction: If a > 0, i’j+ 1/2 = 4 +i/ 2 ; otherwise, 
Vj+i/2 — v j*+\/ 2 - Thus for advection, we need only one of the two values a,l( l r j+i/ 2 - l°r the 
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Euler equations, however, we will need both, and we employ well-known methods for the upwind 
step. 

Equation (1.3) can be integrated by a standard Runge-Kutta method. Here we use the three- 
stage scheme of Shu and Oslier [15]. With v representing { v 3 }, denote by L(v) the spatial operator 

L( v )j ~ —(Vj+i /2 “ v j- 1 / 2 )* (1.4) 

Then this scheme is given by 

= v n 

w?! 1 ) = u;(°) -f aL(w^) 

u;( 2 ) = |u> (0) + + cjL(w {}) )) (1.5) 

wW — -f |(u;( 2 ) + aL(wW)) 

U n +1 — -y;( 3 ) 

where a = ar/h is the CFL number. 

Observe that Runge-Kutta schemes like (1.5) are made up of repeated applications of a single 
stage scheme given by w^ k) + aL(w ik) ). Ar = 0,1, and 2. Moreover, each stage is an explicit Euler 
scheme, e.g., 

W j ' = V J - V J + 1 /2 - Vj- 1 / 2 )- (1 ■<>) 

Therefore, we first design a monotonicity-preserving scheme for (1.6) and then extend it to the full 
scheme (1.5). 


2. Reconstruction 

Without loss of generality, we discuss the reconstruction only for v^ +x/v i.e., we assume a > 0. 
The reconstruction is carried out in two steps. In the first step an accurate and stable formula is 
used to compute the interface value which is called the original value. In the second step, this value 
is then modified or limited appropriately' to achieve a monotonicity-preserving scheme. 

A straightforward choice for vf +l/2 using the five cell averages ..., v J+2 (the same stencil 
as the third-order ENO scheme) is 

v j+ 1/2 = (2t’j -2 ~ + 47vj + 27w J+ i - 3wy +2 )/60. (2.1) 

Other choices include a low phase error fourth-order formula [8] 

vf+ 1/2 = (9»>-2 - 56i?j_i + 194 Vj + 104® j+ , - llv j+2 )/240, (2.2) 
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or a fifth-order accurate implicit formula given by 

(3«f_i/2 + rn ’j+i/2 + t7 f+3/2)/ 10 = + 19i ’j + 10u J+ i)/30. (2.3) 

The implicit formula has the advantage of low dispersive and dissipative errors; its disadvantage is 
that the tridiagonal matrix inversion costs more. 

2.1 Monotonicity constraint 

We derive constraints for the interface value so that monotonicity is preserved by (1.6). First, 
we need a few definitions. Let the median of three numbers be the number that lies between the 
other two. Let minmod (x, y) be the median of x, y , and 0. Equivalently, 

minmod (x, y) = \ [sgn(x) + sgn(j/)] min (|ar|, |t/|). (2.4) 

Conversely, the median function can be expressed in terms of minmod, 

median (x, y, z) — x + minmod (y — x, z — x). (2.5) 

The minmod function can be extended to any number of arguments. For k arguments, minmod (. . 
returns the smallest argument if all arguments are positive, the largest if all are negative, and zero 
otherwise. This function can be coded as 

minmod (si,.. .,-//) = smin(|zi|. \z k \) (2.6a) 

where 

5= i(sgn( 2 i) + sgn( 2 2 )) 5 (sgn(si) + 8gn(5 3 ))...^(sgn(c 1 ) + sgn(rfc))|- (2-6*) 

Also denote by I [z \,..., z^] the interval [min(^i,..., z^), max( Zi ,..., Zk)]- 

We can now derive a simple condition that preserves monotonicitv. At interface j — 1/2, suppose 
the value v f_i /2 hes between Vj-\ and v 3 : 

v j- 1/2 £ I [ v j-i ’ v j \■ (^*0 

Next, for the interface j + 1/2, denote 

v VL = Vj + a( v 3 - Uj_i) (2.8) 
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x j -1 x j x j +1 x 

Figure 2.1: Monotonicity-preserving constraint (2.7) and (2.9). 


where UL stands for upper limit, and a > 2 (more on a momentarily). Suppose the value r l \ . ,,, 
lies between vj and 

4+1/2 € / [v>, l ' UL }- (2.9) 

Then, after one stage via (1.6), the solution lies between iy ] and Vj provided that the time 
step satisfies the condition 

° < 1/(1 + Or). (2.10) 

Indeed, for increasing data, (2.7) and (2.9) imply that the steepest slope v UL - Vj- X satisfies 
v l 1 ~ Vj-i < (q + l)(vj - Vj_i ); therefore, (2.10) implies Vj_ j < < Vj. See Fig. 2.1. 

Note that for parabolic reconstruction schemes, a is typically 2. For Runge-Kutta time stepping, 
we find that a = 4 works well, while o = 2 tends to cause stair-casing. With a — 4, (2.10) leads to 
a CFL number o < 0.2 but, in practice, cr = 0.4 still yields nonoscillatory results. 

Next, assume that (2.7) and (2.9) hold for all j, Expression (2.7) with index j replaced by j -f 1 
takes the form 

V j + 1/2 ^ I[vjjVj+x]. ( 2 . 11 ) 

The above and (2.9) result in the condition that ^ GS * n intersection of the two intervals 

I[ v j' v j+ i] an< J I[ v j,v UL ]. One end of this intersection is Vj. The other is the median of v 3 , , 

and v l L , and is denoted by v MI . Using vj as the pivot, v MP can be expressed by using the minmod 
function, 

v Mr = Vj + ni in mod [v J+1 - Vj,a(v } - v J _ i )]. (2.12) 
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Thus, (2.9) and (2.11) imply 


tff+1/3 € I[Vj,V MP ]. 


(2.13) 


The simplest way to satisfy this constraint is to replace the original Kv the median of vj* +l f 2 , 

Vj, and v MP : 

^yVi/2 median (vj‘ +l/2 ,v j ,v MP ). (2.14) 

Expression (2.14) preserves monotonicity in the following sense: under the CFL restriction 
(2.10) if the data {vj} are monotone, then after one stage, are also monotone. This fact 

follows because lies between Vj- 1 and vj for all j . 

The monotonicity-preserving property extends easily to the full scheme (1.5). Indeed, given 
monotone data {^}, we have just shown that are monotone provided the interface values 

are given by (2.14) and the above CFL restriction is satisfied. Since are monotone, the 

quantities { u + aL(vM^)} are also monotone because they result from applying a single stage 
scheme to Next, for each j, w j 2 ' is a combination of Vj and + a( L( w^))j with positive 

weights independent of j. Therefore, {wj 2 '} are monotone as well. Repeating this argument, it 
follows that {u^} are also monotone. Thus, if the data {nj} are monotone, and the interface 
values are obtained by (2.14), then the cell averages at, the next time level are monotone 

under the CFL restriction (2.10). 

The drawback of (2.14) is that near an extremum, it causes accuracy to degenerate to first- 
order. Figures 2.2(a) and 2.2(b) show the loss of accuracy caused by the constraints (2.11) and 
(2.9), respectively. Note that the data are on a parabola. 

2.2 Accuracy-preserving constraint 

To avoid the loss of accuracy, we enlarge the intervals in (2.11) and (2.9) in such a way that 
these intervals remain the same for monotone data but, near an extremum, these intervals are 
larger, and both contain the original v^ l / 2 . 

First, the interval in (2.11) is enlarged by adjoining the value v MD defined below (MD stands 
for median). At interface j + 1/2, let v FL and v FR be the values extrapolated linearly from the left 
and right, respectively, 

V FL = Vj + Vj - Vj- 1 ), v FR = Vj+i + i(w J+1 - v J+2 ). (2.15) 
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Figure 2.2: Loss of accuracy near extrema: (a ) by (2.11) and (b) by (2.9). 


With 

v ^ — \i v 3 + Vj+x ), (2.16) 

where AV stands for average, set 

t ,A/D _ mec Ji an ( v AV ^ v 1 ' v FR ). (2.17) 

Constraint (2.11) is relaxed to 

v j+ 1/2 £ ? [ v ji v j+u (2.18) 

One can verify that if the four pieces of data {vj -\, vj, v J+ \, v j+2] are monotone, then, at the 
interface j + 1/2, lies between and Vj+i, and the above constraint reduces to (2.11). Near 
an extremum as in the case of Fig. 2.3 (a), however, v AfD lies outside I [vj, ?; J + 1 ] and provides room 
so that the interval I[vj,Vj+ contains the original 

The argument (2.15)-(2.18) conveys the idea. For the purpose of coding, it is more efficient if 
we employ the second differences. Set 


dj - v J _ 1 + Vj+i - 2 vj. 


(2.19) 


and 


^j +)/2 = minmod (dj, cl j+l ), (2.20) 

where MM stands for minmod. Then, since v Fh = v AV - \dj and similarly for v FR , it follows that 

v ^ = v AV_ ¥ ^ /r ( 2 . 21 ) 
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(a) 



(b) 


Figure 2.3: Enlarged monotonicity interval (a) by (2.18) and (b) by (2.23). 


Next, the interval in (2.9) is enlarged by adjoining the value v LC defined below (LC stands for 
large curvature). Consider the parabola p determined by the cell averages Vj-\, vj, and the second 
difference d (a quantity similar to dj). A straightforward calculation gives the value at x J + x / 2 , 

P( x j+ 1 / 2 ) - v j + \( v j ~ Vj -\) + § d. 

The parabola with d — 2 d^Jfj 2 gives 

V LC = Vj + \( Vj - Vj ( 2 - 22 ) 


Constraint (2.9) is relaxed to 


V J+l/2 


e /[»>, 


ML „LC l 


(2.23) 


Using the fact that a > 2, one can verify that if the four pieces of data Vj-i, vj* 

are monotone, then, at the interface j + 1/2, v LC lies between Vj and v L [ \ The above constraint 
therefore reduces to (2.9). Near an extremum as in the case of Fig. 2.3 (b), however, v LC lies outside 
I\vj, v ^L] and provides room so that the interval I [v 3 , v^ ( ] contains the original 

In practice, we reduce the amount of room so that near a non-monotone discontinuity (such 
as a sawtooth profile), constraints (2.18) and (2.23) reduce to (2.11) and (2.9) respectively. This 
reduction can be accomplished by replacing by 


dj+-\/ 2 — minmod (4 dj — dj+i,4dj+i - dj,dj,d/+ 1 )* (2.24) 

To clarify the role of d M ‘ l , assume that d 3 and d J+ x are of t he same sign. Then if dj/dj+x < 1/4 
or dj/dj+i > 4, the above minmod of four arguments returns 0. Loosely put, when the second 


9 







differences change substantially, then d^ /2 = 0, and the extended intervals reduce to the simple 
ones in (2.9) and (2.11). We thus replace expressions (2.21) and (2.22) for each interface j + 1/2 
by 

v md = v av_ ¥ ^ (2.25) 

l ' LC = v j + \( v j - ®>-i) + i d j r -\j 2 ■ (2.26) 

The intersection [v m,n , n max ] of the two intervals I[vj,v J+l ,v MD ) and / [ it ,-, v ul , v lc ] can be calcu- 
lated by 

v mm - inax[ min(wj, v j+l , v MD ), min(t>j, v UL , v LC ) ], (2.27a) 

u max = min[ max( vj, v J+i , v MD ), max( v r v VL , v LC ) ]. (2.276) 

Finally, to ensure that <f +1/2 lies in [u min , u max ], we replace ' J ]+l/2 by the median of i’j' +1/2 , 
t- min , and i’ max : 

v j+i/2 ~ median (i?^ +1/2 , u m,n , t» max ). ( 2 . 28 ) 

The above limiting procedure preserves monotonicity and accuracy. In addition, for most ceils 
in smooth regions, the original sa ^ s ^ GS constraint (2.13) a priori. In this case, the limiting 

procedure (2.24)-(2.28) does not alter the original As a result, we can use (2.13) to detect 

such cells and bypass the limiting procedure altogether. The condition that the original : / 2 

lies in the interval /[ vj,v MP ] is equivalent to (ej +1/2 - ~ v XfP ) < 0. In practice, this 

condition is coded with a tolerance value of c = 10~ 10 : 

( v j+ 1/2 ~ fj)(*»/+ 1/2 “ vMP ) < <■ (2-29) 

We summarize the computation of the interface value below. 

Algorithm for the interface value. Suppose the cell averages {r y } are given, and a > 0. For 
each interface j + 1/2, calculate the original value v p +l/2 from (2.1), and v MP from (2.12). If (2.29) 
holds , then a/lc ^ wc move on t° the next interface. Otherwise, calculate d 3 _ j, d ; , 

d J+1 from (2.19), df + \ /2 and df_\ /2 from (2.24), v UL from (2.8), v AV from (2.16), v MD and v LC 
from (2.25) and (2.26), and v min , t’ max from (2.27). Finally, calculate t:| +1/2 from (2.28), and the 
interface value is then v J+l / 2 = v ,J + l / 2 - 
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2.3 Remarks 

Note that constraint (2.13) on the interface value is a sufficient condition for monotone data 
to remain monotone under Runge-Kutta time stepping. It may be viewed as the analogue of Van 
Leer’s constraint [19] which provides the same type of condition for monotonicity under exact time 
evolution. Also note that the geometric framework, the use of the median function, and v MD were 
introduced by Huynh in [5]. 

For advection with a < 0, the interface value vf+xj 2 1S obtained by reflecting the above expres¬ 
sions about Zj+i/ 2 - To be specific, the reconstruction algorithm is the same with {vj_ 2 , Uj_i, v,, Vj+i, 
Vj+ 2 } replaced by {i>j +3 , Vj +2 ,Vj+i } respectively. Next, the stencil for computing both 

v Hi /2 anc ^ v f+\/2 cons i sts the Slx points ... 15 + 3 }. Therefore, we could define both Vj+i/ 2 

and v j i + i/ 2 ^ the quintic fit of all six cell averages without enlarging the stencil (in the case of 
Euler equations). The corresponding limited scheme, however, is prone to stair-casing. 

Higher-order schemes can be derived using larger stencils. With the same stencil as the /nth- 
order ENO scheme, a (2m - l)th-order scheme can be obtained. For example, for m = 4. we have 


the seven-point formula 


Vj + 1/2 = (—3t>j_3 + 25^-2 - lOli’j-i + 319 Vj 
-f214i?j_i — 38t> ; + 2 + 4 u >+3 )/420, 


and, for m = 5, the nine-point formula 


t’| +1/2 = (4 Vj-4 - 41t?j_ 3 + 199^-2 - 6411V-1 + 1879 w ; - 
+ 1375vj + i - 305 w j+ 2 + 55^+3 - 5t.’ J+ 4 )/2520. 


(2.30a) 


(2.306) 


The same limiting can be employed for these original interface values. The resulting schemes achieve 
high spatial accuracy but remain third-order in time. For the fourth- and fifth-order Runge-Kutta 
methods, in order to preserve monotonicity we need the calculation of the time-reversed operator 
L [15], which is beyond the scope of this paper. 

The above reconstruction depends continuously on the data in the sense that a small change in 
the data causes a small change in the interface value. This property is shared by WENO (but not 


by ENO) reconstruction. 


3. Extensions 


In this section, we describe the extensions of the above schemes to the fuller equations. While 
these extensions are standard, the monotonicity-preserving property may not hold because the equa- 
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tions are nonlinear. Nevertheless, the numerical solutions obtained below are generally nonoscilla- 
tory. 

3.1 Euler system in one dimension 

The Euler equations of gas dynamics for a polytropic gas can be written as 

u,+f(u) x = 0 (3.1) 

where 

u = (p, pu, E) T , 

f(u) = uu + ( 0 ,p, up) T , ( 3 . 2 ) 

P =(7-1 )(E-\pu 2 ). 

Here, T represents the transpose; />, u, p, and E are the density, velocity, pressure, and total energy 
respectively; and 7 = 1.4, is the ratio of specific heats. The speed of sound c is given by (7 p/p) 1 / 2 . 

The eigenvalues of the Jacobian matrix A(u) = di/du are u - c, u and u + c. The matrices of 
left and right eigenvectors of A are needed in the reconstruction. These are given bv ([ 3 ]) 

/ b-ij'l + u/ 2 c —b\u/'2 — 1 / 2 c 6,/2 \ 

L= 1 - b 2 b\ u -61 ( 3 . 3 ) 

\ 62/2 - u/‘2c —byu/2 + l/ 2 c bi/2 J 

and 

1 1 1 \ 

R= u-c u u + c (3.1) 

y H - uc i u 2 H + tic j 

where bi - (7 - l)/c 2 and b 2 - u 2 l>i/ 2 , H = c 2 /(7 - 1) + \u 2 . 

Integrating (3.1) over the cell [xj_ 1/2 , ^+i/ 2 ] yields 

f/u, 1 r i 

“dT + h l f ( u (*i+i/ 2 . 0 ) - f ( u (*j-i/ 2 ./))J = 0 , (3.5) 

where Uj(t) are the cell averages. The first step in calculating f J+1/2 ^ f(u(x J+1/2 , < n )) is to 
reconstruct u on both sides of the interface Xj+i/ 2 - 

It is well known that the reconstruction is best carried out in local characteristic variables [ 3 ]. If 
{vj} are the approximations to the cell averages at time level 11 . these local characteristic variables 
for the cell [zj_i/ 2 , Xj + 1 / 2 ] are given by 

w k = L(Vj)\ J+k , for k = -2,2. (3.0) 
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The scalar reconstruction algorithm is now applied to obtain point values w* 1/2 and wat .r J+1 y 2 
and .Tj_i/ 2 , respectively. The corresponding conservative variables v^ +1/ , 2 and v^_^ 2 are obtained 
by the inverse of (3.6) 

V j + l/2 = R ( v r) w f/2’ v f—1/2 = R ( v j) w -i/2- ( 3 - 7 ) 

At each interface j + 1/2, the two values v^ +1 ^ 2 and 'vf +l j 2 are used to calculate f J+1 / 2 v ia Roe’s 
flux-difference splitting [12]. This splitting is implemented here with Huynh’s entropy fix [6]. 

Equation(3.5) is then integrated by the Runge-Kutta scheme (1.5). The time step is given in 
terms of the CFL number a by 

At = -J--• (3.8) 

Maxjd^l + Cj) 

Note that, the extension described above is standard, but it does not take advantage of the fact 
that our reconstruction algorithm leaves the left and right interface values unchanged in smooth 
regions away from extrema. In these regions, the reconstruction applied to the local characteristic 
variables yields a result identical to formula (2.1) applied to the conserved variables {v ; }. Thus, the 
expense of characteristic decomposition may be avoided for such regions if they could be detected 
in a simple manner as in [6], However, we do not pursue this approach here. 

3.2 Euler system in two dimensions 

An immediate extension of the above scheme to multi-dimensions can be accomplished in the 
same manner as the finite difference ENO schemes of Shu and Osher [15] [16]. The idea is to 
avoid calculating the mixed derivatives of the reconstruction from cell averages by applying the 
reconstruction directly on point values of the fluxes. The reconstruction then reduces to two one¬ 
dimensional reconstructions along the coordinate lines. The same Runge-Kutta scheme (1.5) is used 
to integrate the equations in time. Here we have chosen the Lax-Friedrichs version (ENO-LLF) in 
our numerical experiments. The 2D extension of the schemes derived here are obtained by substi¬ 
tuting our algorithms for reconstruction in place of the scalar one-dimensional ENO reconstruction 
therein. Coding aspects of these schemes can be found in [17]. 

4. Numerical experiments 

For simplicity, we present numerical results only for the scheme combining the quartic fit (2.1) 
and the accuracy-preserving constraint (2.2X). We refer to this scheme as the MP5 scheme (MP 



for monotonicity preserving). Some comparisons with EN03 and WEN05 schemes are provided. 
The three schemes MP5, WEN05 and EN03 have the same stencil, and the first two are spatially 
fifth-order accurate. Listings of the these three reconstruction procedures in Fortran are given in 
Appendix A. Also note that we employ only uniform meshes and, unless otherwise stated, the CEL 
number is 0.4. 

All computations are carried out on a 100 MHz R4000 SGI Indigo Workstation, with c = 10~ 10 
and a = 4. In all cases, the compiler options -r8 -03 were used. We have observed that computing 
times vary widely depending on the hardware and compiler options used. Therefore, computing 
times are to be viewed only as an approximate measure of the efficiency of the various schemes. 
The computing time for the scheme with constant reconstruction and the three-stage Runge-Kutta 
time stepping is also provided as a reference. Since the reconstruction is trivial for this scheme, 
this computing time reflects the cost of all other calculations except reconstruction. 

4.1 Advection of a smooth function 

We solve (1.1) with u(:r,0) = sin(trx) 4 with periodic boundaries. We are particularly interested 
in the behavior of the errors of the cell averages under mesh refinement. Since the function is 
smooth, the most accurate and efficient scheme with the given stencil of five cell averages is the 
unlimited scheme (2.1). We compare the results of the WEN05 and MP5 schemes to this unlimited 
scheme for a - 0.05 in Table 2 (a) and a = 0.4 in Table 2 (b). The results from EX03 are less 
accurate and are not shown. 


Note that the errors obtained by the unlimited scheme and those by the MP5 scheme are 
essentially identical. This confirms that the limiting procedure leaves the quartic fit unchanged 
at. smooth extrema. At low CFL numbers, the MP5 scheme approaches the theoretical order of 
accuracy of five as can be seen in fable 2(a). In both cases, the MP5 scheme compares favorably 
with the WEN05 scheme in both accuracy and efficiency. 
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Table - 2: 

all grids. 


Advection of sin(7rx) 4 by several schemes, t = 2, Ax = 2 /N, CPU time quoted is for 


Table 2(a): At/Ax = 0.05, CPU time for constant reconstruction: 4.30 sec. 


Scheme 



Table 2(b): 


WEN 05 


Too error 


2.3 

3.4 

3.5 
3.4 
1.1 
1.17(1 
1.40(-2 
5.05(-4 
1.63(-5 
5.25(-7 
1.17(1 
1.40(-2 
5.05(-4 
1.63(-5 
5.2 


order 


- 

1.07( 

2.79 

1.73( 

3.29 

1.75( 

3.35 

8.88( 

4.90 

2.54( 

- 

8.05( 

3.06 

8.14( 

4.80 

3 .01 ( 

4.96 

9.74( 

4.95 

3.14( 

- 

00 

o 

3.06 

8.14( 

4.80 

3 .01 ( 

4.96 

9.74( 

4.95 

3 . 14 ( 



L\ order CPU time - sec. 


22.33 


0 .4, CPU time for constant reconstruction: l.Oo sec. 


2.39(1) 
3.74(-2) 
3.26(-3) 
3.00(-4) 
1.25(-5) 
> 1 ) 
- 2 ) 
-3) 



Loc> order 

L\ error 

L\ order 

- 

1-07(1) 

- 

2.68 

1.87(-2) 

2.52 

3.52 

1.79(-3) 

3.39 

3.44 

1.11 (-4) 

4.01 

4.58 

6.17(-6) 

4.17 

- 

8.01(-2) 

- 

2.77 

1.03( 2) 

2.96 

4.01 

6.15(-4) 

4.06 

3.54 

5.05(-5) 

3.61 

3.19 

5.42(-6) 

3.22 

- 

8.01(-2) 

- 

2.77 

1.03(-2) 

2.96 

4.01 

6.17(-4) 

4.06 

3.54 

5.04(-5) 

3.61 

3.19 

5.42(-6) 

3.22 















4.2 Advection of a piecewise continuous function 

Next, the initial condition is given by 

n 0 (yr) = exp(-log(2)(x + 0.7) 2 /0.0009) if - 0.8 < x < -0.6, 

uq(x) = 1 if - 0.4 < x < -0.2, 

a 0 (x) = 1 - |10(x -0.1)| if 0 < x < 0.2, (4-1) 

U 0 (x) = [1 - 100(1 - 0.5) 2 ] 1/2 if 0.4 <x <0.6, 

Uo(x) = 0 otherwise. 

This initial condition includes a Gaussian wave, a square wave, a triangular wave, and a semi¬ 
ellipse. We use 200 cells with a - 0.4. The solutions at t = 2 (after one period or 200 cells) and 
t = 20 (ten periods) are shown in Figs. 4.1 and 4.2 respectively. The solid line represents the exact 
solution. Also shown are the computing times of the various schemes. Again, note that the MP5 
solutions compare well with those by EN03 and WEN05 schemes. 

Resolution at discontinuities can be enhanced by using steepening techniques as in [4], [21]. and 
[6]. These techniques are expensive and, while they are effective in one dimension, it is still not 
clear how well they perform in multi-dimensions. Here, we will limit, our study to the base schemes 
only. 

4.3 Euler system in one dimension 

In the following three problems, the CFL number is 0.4, and the spatial domain is [-1. 1], 
For the initial conditions, unless otherwise stated, the subscript /, denotes -1 < x < 0, and R. 
0 < x < 1. The final time is denoted by </, and the total number of cells, by N. 

1. Sod’s problem [13] 

(Pi,,ul,pl) =(1,0,1), 

(Ph,ur,Pr) = (0.125,0,0.1). 
tf = 0.4, 

N = 100 . 

Since this problem starts from a singularity, smaller time steps are used initially as described in 
[6]. The density field from the MP5 scheme is shown in Fig. 4.3. Note that the contact discontinuity 
and the shock are resolved with high resolution. 
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Figure 4.1: Advection over a period by (a) MP5, (b) WKN05, (c) EN03, (d) Unlimited scheme 
At/Ax = 0.4, Ax = 2/200, t = 2, CPU time for constant reconstruction = 0.03 sec. 
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Figure 4.2: Advection over 10 periods by (a) MP5, (b) WEN05, (c) EN03, (d) Unlimited scheme 
At/Ax = 0.4, Ax = 2/200, t = 20, CPU time for constant reconstruction = 2.09 sec. 
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2. Lax’s problem [10] 


(PL,t£ L ,PL) = (0-445, 0.698,3.528), 

(/>/?, ur,p/?) = (0.5,0,0.571), 

t f = 0.32, 

AT = 100. 

The density field from the MP5 scheme is shown in Fig. 4.4. Again, the contact discontinuity 
and the shock are well-resolved. 

3. Shu’s problem [14] 

In this problem, a moving shock wave interacts with a density disturbance and generates a flow 

field with both smooth structure and discontinuities. Here, L stands for — 1 < x < —0.8, and R , 

— 0.8 < x < 1. The initial conditions, final time, and number of mesh points are 

(PLiULiPl) = (3.857143,2.629369,10.3333), 

(/>/*, = (1 + 0.2sin(57rx),0, 1), 

tj — 0.36, 

N = 300. 

Since the exact solution is not known, the solution by EN03 with 800 cells is used in its place. 

The results of MP5 and WEN05 are shown in Fig. 4.5. The MP5 scheme captures the shock 
with high resolution and resolves all local extrema accurately. 

4.1 Euler system in two dimensions 

We present results for two well known problems. 

1. Oblique shock reflection [1] 

The domain [0,4] x [0,1] is covered by a uniform mesh of 60 x 20 cells. The boundary conditions 
are: at the bottom, solid boundary; at the right, supersonic outflow; at the left, the conditions are 
fixed with 

( p,u,v,p ) = (1,2.9,0,1/7); 

and at the top, 

(/>, !/, v,p) = (1.69997,2.61934, -0.50632,1.52819). 

Under these conditions, an oblique shock forms from the top left corner and is reflected by the 
bottom boundary. Initially, flow conditions at the left boundary are set throughout the whole 
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domain. After 10,000 iterations, the solution essentially reaches steady state. The residual drops 
roughly 2 orders of magnitude for MP5, WEN05, and EN03, while for the minmod and first-order 
upwind schemes the residual drops to machine zero. (Note that the minmod scheme is defined by 
w i /2 = w 0 + 5 minmod (wj - w 0 , w 0 - w_x), where w is the characteristic flux in this case.) 

The pressure along the line y = 0.55 (j=ll) is shown in Fig. 4.6. Concerning accuracy, it 
can be seen that the higher-order schemes have small oscillations about the exact solution. These 
oscillations are reduced on finer grids for all three schemes. Note that the MP5 scheme yields a 
highly accurate solution. 

The computing times of the various schemes are also shown in Fig. 4.6. The first-order and 
minmod schemes are coded here with the local characteristic decompositions over the full five-point 
stencil. This first-order scheme represents the most efficient reconstruction in this framework, and 
its CPU time reflects the overhead of the characteristic decomposition. It can be seen from Fig. 4.6 
that the overhead is more than half of the total computing time. In other words, unlike the case of 
advection, the computing time of the reconstruction step for the Euler equations is less than one 
third of the total time. 

2. Double Mach reflection [20] 

The computational domain is [0,4] X [0,1], The reflecting wall is from (1/6, 0) to (4, 0). Initially, 
a Mach 10 shock is incident on this wall at (1/6,0) making an angle of sixty degrees with the x-axis. 
To the right of the shock is undisturbed fluid of uniform pressure 1 and density 1.4. To the left of 
the shock, the conditions are 

(p. u , v,p) = (8.0, 7.1447. -4.125,116.5). 

As the shock reflects off the wall, a diffraction pattern is formed. The final time is tr = 0.2. A 
detailed description of the problem and various solutions can be found in [20]. 

The boundary conditions are: at the bottom, from (0,0) to (1/6,0), linear extrapolation; from 
(1/6,0) to (4,0), solid boundary; at the right, linear extrapolation; at the left, supersonic inflow; 
at the top, time-dependent conditions determined by the exact motion of the Mach 10 shock. 

The MP5 and WEN05 solutions, obtained using a 240 x 60 mesh, are shown in Fig. 4.7. It can 
be seen that both schemes capture all the significant features of t he solution such as the two Mach 
stems and the wall jet. 
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5. Conclusions 


A new class of high-order schemes for the numerical solution of hyperbolic conservation laws was 
introduced. The key feature of these schemes is the reconstruction procedure which combines an 
accurate interface formula with a monotonicity-preserving constraint. The constraint is designed 
to preserve accuracy and to work well with Runge-Kutta time stepping. It is shown that, for 
advection, if the data are monotone, then the solution is also monotone under a time step restriction. 
Numerical experiments confirm that the resulting scheme is accurate in smooth regions, resolves 
discontinuities with high resolution, and is also efficient. The new scheme compares favorably with 
state-of-the-art schemes such as EN03 and WEN05. 

Appendix-A 


In this appendix, we give the listings of the three higher-order reconstruction algorithms in For¬ 
tran. V(J) are the cell averages vj and VL(J) are the computed interface values ^ +1 / 2 * DMM(X,Y) 
is the minmod function of two arguments while DM4(W,X,Y,Z) is the minmod function of four 
arguments. 


c-MPS RECONSTRUCTION 

c-DMM(X.Y) = 0.5*(3IGN( 1 ,X ) + SIGN(1..Y))« 

Ac MIN ( A B S ( X ), A B S (Y)) 

c-D M 4 (\V, X, Y, Z ) = 0 125*(SIGN(1.,W) + SIGN(1..X))* 

Ac ABS( ( SIGN( 1W ) + 5IGN( 1 -,Y ))* 

Ac (SIGN(1., W ) + SIGN( 1 ,Z)) ) 

A' * MIN ( A B S ( W ), A B S ( X ), A H S (Y), A B S (Z)) 

c 

B 1 = 0.016666666667 
B2 = 1.333333.333333 

ALPHA - 4 
EP5M = l.E-10 
c 

VOR = Bl*(2.*V(J-2)-13.*V(J-l) 

A e ■+ 4 7, * V ( J ) + 27.*V(J + 1) 

A: - 3.*V(J + 2) ) 

VMPs V(J) + DMM(V(J+1)-V(J),ALPHA*(V(J)-V(J-1))) 
IF((VOR*V(J))*(VOR-VMP).LE. EPSM) THEN 
VL( J) = VOR 
ELSE 
ELSE 
C 

DJM1 = V( J-2)-2*V( J-l) + V(J ) 

DJ = V(J-1)-2.*V(J )+ V(J + 1) 

DJP1 = V(J )-2."V(J+l)+ V(J + 2) 

DM4 J P H = DM4(4.*DJ-DJP1,4.*DJP1-DJ,DJ,DJP1) 

DM4 JMH = DM4(4.*DJ-DJM1,4*DJM1-DJ,DJ,DJM1) 

VUL = V(J) + ALPHA “(V(J)-V(J-l)) 

VAV = 0.5*(V(J) + V{J + 1)) 

VMD = VAV - 0 5*DM4JPH 

VLC = V(J) + 0.5*(V( J)-V(J-l)) + B2*DM4JMH 
VMIN = MAX(MIN(V( J) ( V(J + 1),VMD), 

A: MIN(V(J),VUL,VLC)) 

VM AX = MlN(MAX(V(J),V(J + l) t VMD), 

Ac M AX(V(.1), V U L, V LC )) 

VL( J) = VOR + DMMf VMIN-VOR.VMAX-VOR) 

ENDIF 

C-WEN05 RECONSTRUCTION 
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o o 


EPS W = l.E-6 

B1 = 1.083333333333 

B2 = 0.166666666667 

DJM1 = V(J-2)-2 *V(J-1)+ V(J ) 

EJMl = V(J-2)-4 *V(J-1) + 3.*V(J ) 

DJ = V( J-l)-2 *V( J ) + V(J + 1) 

EJ = V(J-l)- V(J+1) 

DJP1 = V(J )-2.*V(J+l) + V(J + 2) 

EJP1 = 3*V(J)-4 *V(J + l)+V(J+2) 
c 

DISO = B1*DJM1*DJM1 + 0.25*EJM 1*EJM 1 + EPSVV 
DIS1 = B1*DJ*DJ + 0.25*EJ*EJ + EPSW 
DIS2 = Bl*DJPl*DJPl + 0 25*EJPl*EJPl + EPSW 
C 

Q30 = 2*V(J-2)-7.*V(J-l) + ll.*V(J ) 

Q31 = -V(J-1) + 5.*V(J )+ 2.*V( J + l) 

Q32 = 2 *V(J )+5.’*V(J +1) -V(J + 2) 

C 

DOl = DIS0/DIS1 

D02 = DIS0/DIS2 

A1BA0 = 6*D01*D01 

A2BA0 = 3.*D02*D02 

WO = 1/(1. + A1B AO + A2BA0) 

W1 - A1 BAO*WO 
W2 = 1. - WO - W1 

VL(J) = B2*( W0*Q30 + W1*Q31 + W2*Q32 ) 

- EN03 RECONSTRUCTION 

DATA CM( 1,1 ),CM{ l,2),CM(l t 3)/2.,-7. 1 ll./ 

DATA CM(2,1 ),CM(2,2),CM(2,3)/-1 .,5.,2./ 

DATA CM( 3,1), C M{ 3,2), CM( 3,3)/2. ,5 1 / 

B 1 = 0.166666666667 

SP = AB5( V(J + 1) - V( J ) ) 

SM = ABS( V(J ) - V(J-l) ) 

DJ = A BS( V(J + 1) - 2. * V { J ) + V( J-l) ) 

IE(2.*SP GT. SM) THEN 

DJM1 = ABS( V(J ) - 2 *V( J-1) + V(J-2) ) 
IF(DJ.GT.2 *DJMl) THEN 
ID = 1 
ELSE 
ID = 2 
ENDIF 
ELSE 

DJP1 = ABS( V( J + 2)-2.*V(J + l) + V(,J ) ) 
IF(2.*DJP1 GT. DJ) THEN 
ID = 2 
ELSE 
ENDIF 
ENDIF 

VL(J)= { CM(ID,1 )*V( J-3 + ID) + CM{ ID,2 )"V( J-2 + ID ) + 
& CM(ID,3)*V( J-l+ID) )*Bl 
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Density, a = 0.4, 100 points 


Density, a = 0.4, 100 points 




Density, a = 0.4, 300 points 



Figure 4.5: Shu's problem with MP5 and WEN05 schemes. 
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Figure 4.6: Oblique shock reflection problem with various schemes, a — 0.4, 60 X 40 grid 
CPU time in seconds is shown above each plot. 
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Figure 4.7(a): Double Mach reflection problem with MP5, CFL = 0.4, 240 X 60 grid. 
30 density contours from 1.73 to 21. CPU time = 5553 seconds. 

CPU time for constant reconst ruction = 4122 seconds. 
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Figure 4.7(b): Double Mach reflection problem with WEN05, CFL = 0.4, 240 X 60 grid. 
30 density contours from 1.73 to 21. CPU time = 6843 seconds. 

CPU time for constant reconstruction= 4122 seconds. 
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